
#include "MCpolygonizer.h"

MCPolygonizer::MCPolygonizer(const int& nvX, const int& nvY, const int& nvZ) {
   // Define the number of cubes
   nvoxelsX = nvX; nvoxelsY = nvY; nvoxelsZ = nvZ;
}

MCPolygonizer::~MCPolygonizer() {
   
   // clear the data
   map<int,int*>::iterator iter = CrossingCells.begin();
   while (iter != CrossingCells.end()) {
      delete [] (*iter).second;
      ++iter;
   }
   CrossingCells.clear();
}

void MCPolygonizer::setCubeValues(int i, GRIDCELL &gc, const vector<Point3> &grid, const vector<Number> &fvalues) {
   
   int ids[8];
   ids[0] = i;
   ids[1] = getIdX(i);
   ids[2] = getIdY(ids[1]);
   ids[3] = getIdY(i);
   ids[4] = getIdZ(i);
   ids[5] = getIdZ(ids[1]);
   ids[6] = getIdZ(ids[2]);
   ids[7] = getIdZ(ids[3]);
   
   for (int j=0; j<8; j++) {
      const Point3 &p = grid[ids[j]];
      const Number &f = fvalues[ids[j]];
      gc.p[j].x = p[0];
      gc.p[j].y = p[1];
      gc.p[j].z = p[2];
      gc.val[j] = f;
      gc.g[j].x = gc.g[j].y = gc.g[j].z = 0;
   }
}

void MCPolygonizer::polygonize(const vector<Point3> &grid, const vector<Number> &fvalues) {
   
   int isolevel = 0;
   float triangles[18*5];
   int indT[15];
   int* newCell=0;
   GRIDCELL gc;
   
   for (int iz=0; iz<nvoxelsZ-1; ++iz) {
      for (int iy=0; iy<nvoxelsY-1; ++iy) {
         for (int ix=0; ix<nvoxelsX-1; ++ix) {
            int id = getCubeId(ix, iy, iz);
            setCubeValues(id, gc, grid, fvalues);
            
            int ntri=0; // number of triangles
            genTriangle(gc, isolevel, &ntri, triangles, indT);
            if (ntri == 0) continue;
            
            newCell = new int[12];
            for (int e=0; e<12; ++e) newCell[e] = -1;
            for (int t=0; t<ntri; ++t) {
               for (int v=0; v<3; ++v) {
                  int indVert = checkNewVertex(ix, iy, iz, indT[(t*3)+v], newCell);
                  if (indVert < 0) {
                     indVert = vertices.size();
                     Point3 p;
                     p[0] = triangles[(t*18)+(v*3)+0];
                     p[1] = triangles[(t*18)+(v*3)+1];
                     p[2] = triangles[(t*18)+(v*3)+2];
                     vertices.push_back(p);
                  }
                  facets.push_back(indVert);
                  newCell[indT[(t*3)+v]] = indVert;
               }
            }
            CrossingCells.insert(make_pair(id, newCell));
         }
      }
   }
}

int MCPolygonizer::checkNewVertex(int i, int j, int k, int edge, int newCell[12]) {

static int cellNeighboorhood[12][9] = {
   { 0,   0,   -1,        0,   -1,   -1,        0,   -1,    0},   // Edge 0
   { 0,   0,   -1,        1,    0,   -1,        1,    0,    0},   // Edge 1
   { 0,   0,   -1,        0,    1,   -1,        0,    1,    0},   // Edge 2
   { 0,   0,   -1,       -1,    0,   -1,       -1,    0,    0},   // Edge 3
   { 0,   0,    1,        0,   -1,    1,        0,   -1,    0},   // Edge 4
   { 0,   0,    1,        1,    0,    1,        1,    0,    0},   // Edge 5
   { 0,   0,    1,        0,    1,    1,        0,    1,    0},   // Edge 6
   { 0,   0,    1,       -1,    0,    1,       -1,    0,    0},   // Edge 7
   {-1,   0,    0,       -1,   -1,    0,        0,   -1,    0},   // Edge 8
   { 1,   0,    0,        1,   -1,    0,        0,   -1,    0},   // Edge 9
   { 1,   0,    0,        1,    1,    0,        0,    1,    0},   // Edge 10
   {-1,   0,    0,       -1,    1,    0,        0,    1,    0}    // Edge 11
   };

static int edgeNeighboorhood[12][3] = {
   {4,    6,    2},  // Edge 0
   {5,    7,    3},  // Edge 1
   {6,    4,    0},  // Edge 2
   {7,    5,    1},  // Edge 3
   {0,    2,    6},  // Edge 4
   {1,    3,    7},  // Edge 5
   {2,    0,    4},  // Edge 6
   {3,    1,    5},  // Edge 7
   {9,   10,   11},  // Edge 8
   {8,   11,   10},  // Edge 9
   {11,    8,   9},  // Edge 10
   {10,    9,   8}   // Edge 11
   };
   
   // se puede optimizar este codigo en la parte del return, ... solo retornar la variavel
   for (int c = 0 ; c < 3 ; c++) { //em funcao de uma aresta, cada cubo tem 3 cubos vizinhos como maximo
      // este IF filtrara solo los cubos vecinos validos para realizar el test
      if (((i+cellNeighboorhood[edge][(c*3)  ]) >= 0) && ((i+cellNeighboorhood[edge][(c*3)  ]) < nvoxelsX-1) && 
          ((j+cellNeighboorhood[edge][(c*3)+1]) >= 0) && ((j+cellNeighboorhood[edge][(c*3)+1]) < nvoxelsY-1) &&
          ((k+cellNeighboorhood[edge][(c*3)+2]) >= 0) && ((k+cellNeighboorhood[edge][(c*3)+2]) < nvoxelsZ-1)) {
         
         int cell = getCubeId(i + cellNeighboorhood[edge][(c*3)],
                        j + cellNeighboorhood[edge][(c*3)+1],
                        k + cellNeighboorhood[edge][(c*3)+2]);
         
         int e = edgeNeighboorhood[edge][c];
         map<int,int*>::const_iterator iter = CrossingCells.find(cell);
         
          if (iter != CrossingCells.end()) {
            if ((*iter).second[e] != -1) {
               return (*iter).second[e];
            }
            else {
               cout<< "+++++++++++++++++++++++++++++++++++++++++++" << endl;
               cout<< "ERROR : Crossing cell vertex not found !!! " << endl;
               cout<< "Cell ( " << i << " , " << j << " , " << k << " ) " << endl;
               cout<< "edge    = " << edge << endl;
               cout<< "Neighboorhood (" <<cellNeighboorhood[edge][(c*3)  ] << " , " << 
                                          cellNeighboorhood[edge][(c*3)+1] << " , " << 
                                          cellNeighboorhood[edge][(c*3)+2] << " ) " << endl;
                              
               cout<< "neig (" <<i + cellNeighboorhood[edge][(c*3)  ] << " , " << 
                                 j + cellNeighboorhood[edge][(c*3)+1] << " , " << 
                                 k + cellNeighboorhood[edge][(c*3)+2] << " ) " << endl;
               cout<< "cell    = " << cell << endl;
               cout<< "e       = " << e << endl;
               cout<< "c       = " << c << endl;
               cout<< "++++++++++++++++++++++++++++++++++++++++++++" << endl; 
               return -1;
            }
         }
         else {
         }
      }
   }
   return newCell[edge];
}

void MCPolygonizer::export2off(char *filename) {
   
   assert(not facets.empty());
   assert(not vertices.empty());
   
   std::ofstream fileout(filename);
   fileout<<"OFF"<<endl;
   fileout<<vertices.size()<<" "<<facets.size()/3.0<<" "<<"0"<<endl;
   
   for (uint i=0; i<vertices.size(); i++) {
      fileout<<vertices[i][0]<<" "<<vertices[i][1]<<" "<<vertices[i][2]<<endl;
   }
   for (uint i=0; i < facets.size(); i+=3) {
      fileout<<"3"<<" "<<facets[i]<<" "<<facets[i+2]<<" "<<facets[i+1]<<endl;      
   }   
   fileout.close();
}
